# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "patchwork",
"viridis", "RColorBrewer", "here", "ggnewscale", "ggh4x")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
pal_temp <- brewer.pal(n = 10, name = "Paired")[c(2, 6)]
pal_oxy <- brewer.pal(n = 10, name = "Paired")[c(8, 4)]
# Set path
home <- here::here()Plot conditional effects and grid predictions from 02
Intro
Plot conditional effects from predictions in script “02-fit-sdms-conditional.qmd”
Load packages & source functions
Read predictions
preds <- read_csv(paste0(home, "/output/conditional_02_sdms.csv")) %>%
filter(oxy > 3 & oxy < 9) %>%
filter(temp > 2 & temp < 14) %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) %>%
group_by(group) %>%
mutate(est_sc = exp(est)/max(exp(est))) # for plotting all together in heatmapRows: 15000 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): group
dbl (14): temp, oxy, quarter, temp_sc, temp_sq, oxy_sc, year, year_f, quarte...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 5,700 rows (38%), 9,300 rows remaining
filter: removed 2,418 rows (26%), 6,882 rows remaining
mutate: changed 6,882 values (100%) of 'species' (0 new NA)
changed 6,882 values (100%) of 'life_stage' (0 new NA)
group_by: one grouping variable (group)
mutate (grouped): new variable 'est_sc' (double) with 3,658 unique values and 0% NAgrid_preds <- read_csv(paste0(home, "/output/pred_grids_02_sdms.csv")) %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))Rows: 5660568 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): month_year, ices_rect, group
dbl (28): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: changed 5,660,568 values (100%) of 'species' (0 new NA)
changed 5,660,568 values (100%) of 'life_stage' (0 new NA)cogs <- read_csv(paste0(home, "/output/cogs_02_sdms.csv")) %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))Rows: 174 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): group
dbl (9): year, est_x, lwr_x, upr_x, se_x, est_y, lwr_y, upr_y, se_y
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: changed 174 values (100%) of 'species' (0 new NA)
changed 174 values (100%) of 'life_stage' (0 new NA)Calculate weighted quantiles
wq <- grid_preds %>%
group_by(group, quarter) %>%
summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
"10_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.1),
"25_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
"75_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
"90_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.9),
"Env_oxy" = median(oxy),
"Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
"10_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.1),
"25_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
"75_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
"90_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.9),
"Env_temp" = median(temp)) %>%
ungroup() %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) group_by: 2 grouping variables (group, quarter)summarise: now 12 rows and 14 columns, one group variable remaining (group)ungroup: no grouping variablesmutate: changed 12 values (100%) of 'species' (0 new NA) changed 12 values (100%) of 'life_stage' (0 new NA)wq_annual <- grid_preds %>%
group_by(group, year, quarter) %>%
# Could use reframe here...
summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
"1st_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
"3rd_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
"Env_oxy" = median(oxy),
"Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
"1st_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
"3rd_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
"Env_temp" = median(temp)) %>%
ungroup() %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) group_by: 3 grouping variables (group, year, quarter)summarise: now 348 rows and 11 columns, 2 group variables remaining (group, year)ungroup: no grouping variablesmutate: changed 348 values (100%) of 'species' (0 new NA) changed 348 values (100%) of 'life_stage' (0 new NA)Plot conditional effects
# Oxygen
preds_oxy <- preds %>%
filter(temp %in% sort(unique(preds$temp))[0.5*length(unique(preds$temp))])filter (grouped): removed 6,696 rows (97%), 186 rows remainingwi_range_oxy <- preds_oxy %>%
left_join(wq %>%
filter(quarter == 4) %>%
dplyr::select(group, `10_oxy`, `25_oxy`, `75_oxy`, `90_oxy`), by = "group")filter: removed 6 rows (50%), 6 rows remainingleft_join: added 4 columns (10_oxy, 25_oxy, 75_oxy, 90_oxy) > rows only in x 0 > rows only in y ( 0) > matched rows 186 > ===== > rows total 186# TODO: make predictions for the exact value??
wi_median_oxy <- preds_oxy %>%
left_join(wq %>%
filter(quarter == 4) %>%
dplyr::select(group, Median_oxy), by = "group") %>%
filter(oxy > 0.985*Median_oxy & oxy < 1.015*Median_oxy) %>%
as.data.frame()filter: removed 6 rows (50%), 6 rows remainingleft_join: added one column (Median_oxy) > rows only in x 0 > rows only in y ( 0) > matched rows 186 > ===== > rows total 186filter (grouped): removed 180 rows (97%), 6 rows remainingcol_oxy <- "grey30"#viridis(n = 12, option = "mako")[7]
col_temp <- "grey30"#viridis(n = 12, option = "plasma")[8]
p1 <- ggplot() +
ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
# NOTE: 80% Confidence interval
# 25th and 75th percentile
geom_ribbon(data = wi_range_oxy %>% filter(oxy > `25_oxy` & oxy < `75_oxy`), fill = col_oxy,
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_oxy %>% filter(oxy > `25_oxy` & oxy < `75_oxy`), aes(oxy, exp(est)),
color = col_oxy, alpha = 0.5, linewidth = 1.5) +
# 10th and 90th percentile
geom_ribbon(data = wi_range_oxy %>% filter(oxy > `10_oxy` & oxy < `90_oxy`), fill = col_oxy,
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_oxy %>% filter(oxy > `10_oxy` & oxy < `90_oxy`), aes(oxy, exp(est)),
color = col_oxy, alpha = 0.5, linewidth = 1) +
# NOTE: 85% Confidence interval
geom_ribbon(data = preds_oxy, fill = col_oxy,
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = preds_oxy, aes(oxy, exp(est)), alpha = 0.5, linewidth = 0.5, color = col_oxy) +
# Median!
geom_point(data = wi_median_oxy, aes(Median_oxy, exp(est)),
color = col_oxy, alpha = 0.5, size = 2.5) +
labs(x = "Oxygen (ml/L)", y = "Biomass density (kg/km<sup>2</sup>)") +
theme(#legend.position = "bottom",
#legend.spacing.y = unit(0.1, 'cm'),
#legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
axis.title.y = ggtext::element_markdown(),
plot.margin = unit(c(0,0,0,0), "cm")) +
NULLfilter (grouped): removed 141 rows (76%), 45 rows remainingfilter (grouped): removed 141 rows (76%), 45 rows remainingfilter (grouped): removed 101 rows (54%), 85 rows remaining
filter (grouped): removed 101 rows (54%), 85 rows remaining# Temperature
preds_temp <- preds %>%
filter(oxy %in% sort(unique(preds$oxy))[0.5*length(unique(preds$oxy))])filter (grouped): removed 6,660 rows (97%), 222 rows remainingwi_range_temp <- preds_temp %>%
left_join(wq %>%
filter(quarter == 4) %>%
dplyr::select(group, `10_temp`, `25_temp`, `75_temp`, `90_temp`), by = "group")filter: removed 6 rows (50%), 6 rows remainingleft_join: added 4 columns (10_temp, 25_temp, 75_temp, 90_temp) > rows only in x 0 > rows only in y ( 0) > matched rows 222 > ===== > rows total 222# TODO: make predictions for the exact value??
wi_median_temp <- preds_temp %>%
left_join(wq %>%
filter(quarter == 4) %>%
dplyr::select(group, Median_temp), by = "group") %>%
filter(temp > 0.985*Median_temp & temp < 1.015*Median_temp) %>%
as.data.frame()filter: removed 6 rows (50%), 6 rows remainingleft_join: added one column (Median_temp) > rows only in x 0 > rows only in y ( 0) > matched rows 222 > ===== > rows total 222filter (grouped): removed 216 rows (97%), 6 rows remainingwi_median_temp temp oxy quarter temp_sc temp_sq oxy_sc year year_f quarter_f
1 7.698927 5.7852 1 0.6728450 0.4527204 -0.1932578 1999 1999 1
2 8.024739 5.7852 1 0.7909722 0.6256370 -0.1927308 1999 1999 1
3 10.305422 5.7852 1 1.6048748 2.5756232 -0.1762996 1999 1999 1
4 9.653798 5.7852 1 1.3694188 1.8753080 -0.1760833 1999 1999 1
5 8.350551 5.7852 1 0.8781839 0.7712069 -0.2453884 1999 1999 1
6 8.350551 5.7852 1 0.8779631 0.7708193 -0.2440774 1999 1999 1
sal_sc depth_sc depth_sq group species life_stage est
1 0 0 0 cod_adult Cod Adult 3.4410899
2 0 0 0 cod_juvenile Cod Juvenile 1.5199580
3 0 0 0 plaice_adult Plaice Adult -1.4813702
4 0 0 0 plaice_juvenile Plaice Juvenile -6.4589174
5 0 0 0 flounder_adult Flounder Adult 4.7846446
6 0 0 0 flounder_juvenile Flounder Juvenile -0.5982824
est_se est_sc Median_temp
1 0.6233005 0.9122991 7.598575
2 0.8211386 0.6403539 8.116226
3 0.8638716 0.7718770 10.382043
4 1.2361272 0.5520840 9.570449
5 0.3197793 0.4830844 8.440493
6 0.5137882 0.4085891 8.330722p2 <- ggplot() +
ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
# NOTE: 80% Confidence interval
# 25th and 75th percentile
geom_ribbon(data = wi_range_temp %>% filter(temp > `25_temp` & temp < `75_temp`), fill = col_temp,
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_temp %>% filter(temp > `25_temp` & temp < `75_temp`), aes(temp, exp(est)),
color = col_temp, alpha = 0.5, linewidth = 1.5) +
# 10th and 90th percentile
geom_ribbon(data = wi_range_temp %>% filter(temp > `10_temp` & temp < `90_temp`), fill = col_temp,
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_temp %>% filter(temp > `10_temp` & temp < `90_temp`), aes(temp, exp(est)),
color = col_temp, alpha = 0.5, linewidth = 1) +
# NOTE: 85% Confidence interval
geom_ribbon(data = preds_temp, fill = col_temp,
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = preds_temp, aes(temp, exp(est)), alpha = 0.5, linewidth = 0.5, color = col_temp) +
# Median!
geom_point(data = wi_median_temp, aes(Median_temp, exp(est)),
color = col_temp, alpha = 0.5, size = 2.5) +
labs(x = "Temperature (°C)", y = "Biomass density (kg/km<sup>2</sup>)") +
theme(#legend.position = "bottom",
#legend.spacing.y = unit(0.1, 'cm'),
#legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
axis.title.y = ggtext::element_markdown(),
plot.margin = unit(c(0,0,0,0), "cm")) +
NULLfilter (grouped): removed 178 rows (80%), 44 rows remainingfilter (grouped): removed 178 rows (80%), 44 rows remainingfilter (grouped): removed 138 rows (62%), 84 rows remaining
filter (grouped): removed 138 rows (62%), 84 rows remaining(p1 / p2) + plot_annotation(tag_levels = "A")
ggsave(paste0(home, "/figures/oxy_temp_conditional.pdf"), width = 17, height = 17, units = "cm")Plot weighted quantiles over time
t <- wq_annual %>%
pivot_longer(c("Median_oxy", `1st_oxy`, `3rd_oxy`, "Env_oxy", "Median_temp", `1st_temp`, `3rd_temp`, "Env_temp")) %>%
mutate(var = ifelse(grepl("oxy", name), "Oxygen", "Temperature"),
type = ifelse(grepl("Env", name), "Environment", "Biomass-weighted")) %>%
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage),
group2 = paste(species, life_stage, sep = " ")) %>%
filter(name %in% c("Median_oxy", "Median_temp", "Env_oxy", "Env_temp"))pivot_longer: reorganized (Median_oxy, 1st_oxy, 3rd_oxy, Env_oxy, Median_temp, …) into (name, value) [was 348x13, now 2784x7]mutate: new variable 'var' (character) with 2 unique values and 0% NA new variable 'type' (character) with 2 unique values and 0% NAmutate: changed 2,784 values (100%) of 'species' (0 new NA) changed 2,784 values (100%) of 'life_stage' (0 new NA) new variable 'group2' (character) with 6 unique values and 0% NAfilter: removed 1,392 rows (50%), 1,392 rows remainingt_env <- t %>%
filter(type == "Environment" & group == "cod_adult") %>%
mutate(group2 = "Environment")filter: removed 1,276 rows (92%), 116 rows remainingmutate: changed 116 values (100%) of 'group2' (0 new NA)t2 <- bind_rows(t_env, t %>% filter(!type == "Environment"))filter: removed 696 rows (50%), 696 rows remaining# Reorder to have Environment at the end
t2$group2 <- factor(t2$group2, levels = rev(unique(t2$group2)))
pal <- brewer.pal(name = "Dark2", n = 8)[3:8]
t2 %>%
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter4")) %>%
ggplot(aes(year, value, color = group2, linetype = group2)) +
geom_line(alpha = 0.8) +
#facet_grid(quarter ~ var) +
ggh4x::facet_grid2(var ~ quarter2, scales = "free") +
labs(y = "Environmental variable", x = "Year", linetype = "", color = "") +
guides(color = guide_legend(nrow = 3,override.aes = list(linetype = c(rep(1, 6), 2),
size = 0.5)),
linetype = "none") +
scale_x_continuous(breaks = c(seq(1993, 2021, by = 5))) +
scale_linetype_manual(values = c(rep(1, 6), 2)) +
scale_color_manual(values = c(rev(pal), "grey40")) +
theme(legend.position = c(0.24, 0.36),
legend.spacing.y = unit(-2, 'cm'),
legend.spacing.x = unit(0, 'cm'),
legend.background = element_rect(fill = NA))mutate: new variable 'quarter2' (character) with 2 unique values and 0% NA
ggsave(paste0(home, "/figures/wm_temp_oxy.pdf"), width = 17, height = 14, units = "cm")Plot centrs of gravity
plot_map +
geom_point(data = cogs, aes(est_x, est_y, color = year), size = 0.5) +
geom_errorbar(data = cogs, aes(y = est_y, xmin = lwr_x, xmax = upr_x, color = year),
width = 0, alpha = 0.3, linewidth = 0.3) +
geom_errorbar(data = cogs, aes(x = est_x, ymin = lwr_y, ymax = upr_y, color = year),
width = 0, alpha = 0.3, linewidth = 0.3) +
geom_path(data = cogs, aes(est_x, est_y, color = year), linewidth = 0.3) +
facet_grid(life_stage ~ species) +
labs(color = "Year") +
theme(legend.position = c(0.05, 0.85),
legend.key.height = unit(0.3, "cm"),
legend.key.width = unit(0.1, "cm"))
ggsave(paste0(home, "/figures/cogs_map.pdf"), width = 17, height = 11, units = "cm")
pal <- brewer.pal(name = "Dark2", n = 8)[3:8]
cogs <- cogs %>%
mutate(group2 = paste(species, life_stage, sep = " "))mutate: new variable 'group2' (character) with 6 unique values and 0% NAp1 <- ggplot(cogs, aes(year, est_x, color = group2, fill = group2)) +
geom_ribbon(aes(ymin = lwr_x, ymax = upr_x), color = NA, alpha = 0.3) +
geom_line() +
scale_color_manual(values = c(pal)) +
scale_fill_manual(values = c(pal)) +
guides(color = "none", fill = "none") +
labs(x = "", y = "UTM X center of gravity", color = "", fill = "")
p2 <- ggplot(cogs, aes(year, est_y, color = group2, fill = group2)) +
geom_ribbon(aes(ymin = lwr_y, ymax = upr_y), color = NA, alpha = 0.3) +
geom_line() +
scale_color_manual(values = c(pal)) +
scale_fill_manual(values = c(pal)) +
guides(color = guide_legend(nrow = 2)) +
labs(x = "", y = "UTM Y center of gravity", color = "", fill = "") +
theme(legend.position = "bottom",
legend.spacing.y = unit(-2, 'cm'),
legend.spacing.x = unit(0, 'cm'),
legend.background = element_rect(fill = NA))
p1 / p2
ggsave(paste0(home, "/figures/supp/cogs.pdf"), width = 17, height = 20, units = "cm")Correlation and slope plot
# Correlation plot
dcor <- t %>% filter(!type == "Environment") %>%
mutate(id = paste(name, year, quarter, sep = "_")) %>%
dplyr::select(group2, year, var, value, quarter, id)filter: removed 696 rows (50%), 696 rows remainingmutate: new variable 'id' (character) with 116 unique values and 0% NAt_env2 <- t_env %>%
mutate(name = ifelse(name == "Env_oxy", "Median_oxy", "Median_temp")) %>%
mutate(id = paste(name, year, quarter, sep = "_")) %>%
rename(env_value = value) %>%
dplyr::select(id, env_value)mutate: changed 116 values (100%) of 'name' (0 new NA)
mutate: new variable 'id' (character) with 116 unique values and 0% NArename: renamed one variable (env_value)dcor2 <- dcor %>% left_join(t_env2, by = "id")left_join: added one column (env_value) > rows only in x 0 > rows only in y ( 0) > matched rows 696 > ===== > rows total 696cors <- plyr::ddply(dcor2, c("group2", "var", "quarter"),
summarise, cor = round(cor(env_value, value), 2)) %>%
arrange(var)summarise: now one row and one column, ungroupedsummarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungroupedhist(cors$cor)
# Plot correlations between environment and weighted
cor_pal <- brewer.pal(n = 8, name = "Dark2")[6:7]
dcor2 %>%
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4")) %>%
ggplot(aes(env_value, value, color = quarter2)) +
ggh4x::facet_grid2(var ~ group2, scales = "free") +
geom_abline(color = "grey30", linetype = 2, linewidth = 0.5) +
geom_point(alpha = 0.7, size = 0.5) +
geom_text(data = cors %>% filter(quarter == 1),
aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 1.5,
inherit.aes = FALSE, fontface = "italic", color = cor_pal[2], size = 2.5) +
geom_text(data = cors %>% filter(quarter == 4),
aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 2.7,
inherit.aes = FALSE, fontface = "italic", color = cor_pal[1], size = 2.5) +
labs(y = "Biomass-weighted", x = "Environment", color = "") +
stat_smooth(method = "lm", se = FALSE, linewidth = 0.5) +
scale_color_manual(values = cor_pal) +
scale_x_continuous(breaks = c(4:8)) +
theme_sleek(base_size = 10) +
theme(aspect.ratio = 1,
legend.key.size = unit(0.2, 'cm'),
legend.text = element_text(size = 6),
legend.position = c(0.11, 0),
legend.spacing.y = unit(-1, 'cm'))mutate: new variable 'quarter2' (character) with 2 unique values and 0% NAfilter: removed 12 rows (50%), 12 rows remaining
filter: removed 12 rows (50%), 12 rows remaining`geom_smooth()` using formula = 'y ~ x'
ggsave(paste0(home, "/figures/supp/wm_cors.pdf"), width = 17, height = 7, units = "cm")`geom_smooth()` using formula = 'y ~ x'# Plot the slopes over time for group and quarter, compare to the environment slope. Boxplots and lines??
s_slopes <- t %>% filter(!type == "Environment") %>%
mutate(id = paste(group, var, quarter, sep = "_")) %>%
ungroup() %>%
dplyr::select(year, id, value)filter: removed 696 rows (50%), 696 rows remaining
mutate: new variable 'id' (character) with 24 unique values and 0% NA
ungroup: no grouping variablest_slopes <- t_env %>%
mutate(id = paste("Environment", group2, var, quarter, sep = "_")) %>%
dplyr::select(year, id, value)mutate: new variable 'id' (character) with 4 unique values and 0% NAslopes <- bind_rows(s_slopes, t_slopes)
# Calculate slopes over time
slopes2 <- slopes %>%
split(.$id) %>%
purrr::map(~lm(value ~ year, data = .x)) %>%
purrr::map_df(broom::tidy, .id = 'Year') %>%
filter(term == "year") %>%
mutate(upr = estimate + 1.96*std.error,
lwr = estimate - 1.96*std.error,
sig = ifelse(estimate > lwr & estimate < upr, "sig.", "not sig.")) %>%
rename(id = Year,
year_slope = estimate) %>%
dplyr::select(id, year_slope, sig, lwr, upr) %>%
separate(id, into = c("species", "life_stage", "variable", "quarter"), remove = FALSE) %>%
#mutate(group3 = paste(str_to_title(species), str_to_title(life_stage), paste0("Q", quarter)))
mutate(group3 = paste(str_to_title(species), str_to_title(life_stage)),
x = "x") %>%
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4"))filter: removed 28 rows (50%), 28 rows remaining
mutate: new variable 'upr' (double) with 28 unique values and 0% NA
new variable 'lwr' (double) with 28 unique values and 0% NA
new variable 'sig' (character) with one unique value and 0% NA
rename: renamed 2 variables (id, year_slope)
mutate: new variable 'group3' (character) with 7 unique values and 0% NA
new variable 'x' (character) with one unique value and 0% NA
mutate: new variable 'quarter2' (character) with 2 unique values and 0% NAhlines <- slopes2 %>%
filter(species == "Environment")filter: removed 24 rows (86%), 4 rows remainingpal <- brewer.pal(name = "Dark2", n = 8)[3:8]
slopes2 %>%
filter(!species == "Environment") %>%
ggplot(aes(x, year_slope, ymin = lwr, ymax = upr, color = group3)) +
geom_rect(data = hlines, aes(xmin = -Inf, xmax = Inf, ymin = lwr, ymax = upr, fill = "Env. slope 95% CI"),
inherit.aes = FALSE, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = 3, color = "tomato2", alpha = 1) +
geom_hline(data = hlines, aes(yintercept = year_slope, linetype = "Env. slope"), alpha = 0.7) +
geom_point(size = 2.5, position = position_dodge(width = 0.4)) +
geom_errorbar(position = position_dodge(width = 0.4), width = 0) +
ggh4x::facet_grid2(variable ~ quarter2, scales = "free") +
scale_fill_manual(values = "grey10") +
scale_linetype_manual(values = 2) +
scale_color_manual(values = pal) +
labs(color = "", y = "Year-slope", linetype = "", fill = "") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
legend.position = "bottom") +
NULLfilter: removed 4 rows (14%), 24 rows remaining
ggsave(paste0(home, "/figures/year_slopes.pdf"), width = 17, height = 17, units = "cm")Plot spatially-varying quarter effects, spatial predictions, and ST random fields
grid_preds <- grid_preds %>%
mutate(group2 = paste(species, life_stage, sep = " ")) mutate: new variable 'group2' (character) with 6 unique values and 0% NAnames(grid_preds) [1] "X" "Y" "year"
[4] "lon" "lat" "depth"
[7] "quarter" "month" "month_year"
[10] "oxy" "temp" "sal"
[13] "sub_div" "ices_rect" "temp_sc"
[16] "temp_sq" "oxy_sc" "oxy_sq"
[19] "sal_sc" "depth_sc" "depth_sq"
[22] "quarter_f" "year_f" "est"
[25] "est_non_rf" "est_rf" "omega_s"
[28] "zeta_s_quarter_f1" "zeta_s_quarter_f4" "epsilon_st"
[31] "group" "species" "life_stage"
[34] "group2" # Plot spatially-varying effect
sv <- grid_preds %>%
#filter(group == pull(grid_preds, group)[1]) %>%
filter(year == 1999) %>%
pivot_longer(c("zeta_s_quarter_f1", "zeta_s_quarter_f4"),
names_to = "field", values_to = "zeta") %>%
mutate(Quarter = ifelse(field == "zeta_s_quarter_f1", "Quarter 1", "Quarter 4"))filter: removed 5,465,376 rows (97%), 195,192 rows remainingpivot_longer: reorganized (zeta_s_quarter_f1, zeta_s_quarter_f4) into (field, zeta) [was 195192x34, now 390384x34]mutate: new variable 'Quarter' (character) with 2 unique values and 0% NAplot_map_fc +
geom_raster(data = sv, aes(X*1000, Y*1000, fill = zeta)) +
facet_grid(Quarter ~ group2) +
scale_fill_gradient2(name = "Zeta") +
theme_sleek(base_size = 8) +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"))Warning: Removed 72 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/sv_effects.pdf"), width = 17, height = 9, units = "cm")Warning: Removed 72 rows containing missing values (`geom_raster()`).for(i in unique(grid_preds$group2)){
dd <- grid_preds %>% filter(group2 == i)
j <- pull(dd, group)[1]
# Quarter 1
print(
plot_map_fc +
geom_raster(data = filter(dd, quarter == 1), aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
title = i, subtitle = "Quarter 1")
)
ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q1_", j, ".pdf")), width = 17, height = 17, units = "cm")
# Quarter 4
print(
plot_map_fc +
geom_raster(data = filter(dd, quarter == 4), aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
title = i, subtitle = "Quarter 4")
)
ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q4_", j, ".pdf")), width = 17, height = 17, units = "cm")
}filter: removed 4,717,140 rows (83%), 943,428 rows remainingfilter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remainingWarning: Removed 87 rows containing missing values (`geom_raster()`).Warning: Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 87 rows containing missing values (`geom_raster()`).
Removed 87 rows containing missing values (`geom_raster()`).